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We present high-accuracy calculations of the density of states using multicanonical methods for 
lattice gauge theory with a compact gauge group U(l) on 4 4 , 6 4 and 8 4 lattices. We show that 
the results are consistent with weak and strong coupling expansions. We present methods based on 
Chebyshev interpolations and Cauchy theorem to find the (Fisher's) zeros of the partition function 
in the complex j3 = l/g 2 plane. The results are consistent with reweighting methods whenever the 
latter are accurate. We discuss the volume dependence of the imaginary part of the Fisher's zeros, 
the width and depth of the plaquette distribution at the value of /3 where the two peaks have equal 
height. We discuss strategies to discriminate between first and second order transitions and explore 
them with data at larger volume but lower statistics. Higher statistics and even larger lattices are 
necessary to draw strong conclusions regarding the order of the transition. 

PACS numbers: 11.15.-q, 11. 15. Ha, 11.15.Me, 12.38.Cy 



I. INTRODUCTION 

The ongoing effort at the Large Hadron Collider has 
trigerred a renewed interest in the phase structure of lat- 
tice gauge theory models that may possibly provide an 
alternative to the Higgs mechanism of the electro-weak 
symmetry breaking. In particular, the location of the 
conformal windows for various families of models have 
stirred intense discussions. Different numerical and an- 
alytical techniques have been applied to QCD-like mod- 
els with large number of fermion flavors (Tj-Q or with 
fermions in higher representations 04HJ ■ An interesting 
attempt to classify possible phases of such models based 
on an effective potential for the Polyakov loop was made 
in See also (T34l6| for recent reviews of results and 
expectations. Another direction where massive vector 
bosons emerge without introducing new fermion species 
but in a model with modified gauge transformations has 
been pursued in [It], HH ■ 

In this context, it is important to understand the crit- 
ical behavior of lattice models from as many consistent 
points of view as possible. Recently, it was proposed 
to consider complex extensions jl9l - l2l| of the framework 
proposed by Tomboulis [22j to explain confinement from 



the point of view of the Renormalization Group (RG) 
approach. A general feature that we observed is that the 
Fisher's zeros, the zeros of the partition function in the 
complex /3 plane [23| . play an important role in the de- 
termination of the global properties of the complex RG 
flows. In the case where a phase transition is present, the 
scaling properties of the zeros p^ - [29j allow us to distin- 
guish between a first and second order phase transition. 

Despite its apparent simplicity, the case of the 4D pure 
gauge compact U(l) model with a Wilson action is not 
completely free of controversy. The presence of a double 
peak for the plaquette distributions near /? ~ 1 suggests 
a first order phase transition. However, if spherical or 
toroidal lattices are considered, the double peak disap- 
pears [25|,[26|. In addition, finite size scaling, at relatively 
small volumes seems consistent with a second order phase 
transition with an exponent v ~ 0.35 — 0.40. On the 
other hand, a possible scenario [30| is that as the volume 
increases, v slowly "rolls" towards the first-order value 
i/ = l/£) = 0.25. In the more recent literature (3lTl33 |. 
the idea that the transition is first order is favored. Using 
finite size scaling, the authors of Ref. [33| estimated the 
critical value /3oo = 1.0111310(62). 

In this article, we introduce new methods to locate the 



2 



Fisher's zeros of the the 4D pure gauge compact U(l) 
model with a Wilson action. We rely on high accuracy 
determinations of the density of states, a quantity defined 
in Sec. HH by multicanonical methods |35M38j presented 
in Sec. IIIII The lattice sizes considered are 4 4 , 6 4 and 
8 4 . The consistency of the results with weak and strong 
coupling expansion is checked in Sec. IIV1 The density of 
states has a convex region which implies a double peak 
plaqucttc distribution near ~ 1. The volume depen- 
dence of the double peak distribution is discussed empir- 
ically in Sec. |Vj In Sec. IVI1 Chebyshev interpolations of 
the logarithm of the density of states and Cauchy the- 
orem are used to find the Fisher's zeros in the complex 
j3 = 1/g 2 plane. For the lowest zeros, it is possible to 
check consistency with reweighting methods within er- 
ror bars estimated from the statistical fluctuations of 20 
independent multicanonical streams. 

In the following, we use very high statistics on rather 
small lattices, because this allow us to explore new anal- 
ysis methods and to test whether they converge faster 
towards the infinite volume limit. We plan to use simi- 
lar methods for SU(2) where the imaginary parts of the 
lowest zeros are larger and reweighting methods become 
less reliable when the volume increases (39| . 

Using high statistics at small volumes (4 4 , 6 4 and 8 4 ), 
we show that the imaginary part of the lowest zero ap- 
pears to decrease like L~ 3m , when the linear size L in- 
creases from 4 to 8. This could be interpreted as signaling 
a second-order phase transition with v = 0.325, a value 
close to the estimates of Refs. [25|, |26[. However, using 
data at larger volumes but with lower statistics, we found 
indications for the "rolling" scenario of [3(J . This is dis- 
cussed in Sec. IVIII where we also consider volume effects 
on the width and depth of the plaquette distribution at 
the value of j3 where the two peaks have equal height. 
Simulations required to provide a clear cut distinction 
between first and second order transitions are discussed 
in the Conclusions. 



II. DENSITY OF STATES IN ABELIAN GAUGE 
THEORY 

In the following, we consider the pure gauge partition 
function 

i j 

with /3 = 1/g 2 and the Wilson action 

S = ^{1 - cos9 p ) . (2) 
p 

We use D dimensional symmetric cubic lattices with 
L D sites and periodic boundary conditions. The number 
of plaquettes is denoted M p = L D D(D - l)/2. We 
define the average action: 

P = {S/Np) = -d(ln Z/N p )/df3 . (3) 



Inserting 1 as the integral of the delta function over S 
in Z, we can write 

r 2Af v 

Z = dS n(S) e~P s , (4) 

Jo 

with the density of states defined as 

n(S) - II / dUlS ( S ~ X^ 1 _ cos6 p^ ■ ( 5 ) 
i J p 

Furthermore, we introduce the notation s for S/M p and 
we define the entropy density /(s) via the equation 

n(S) = e^p/W^p) . (6) 

A more general discussion for spin models (40| or gauge 
models [4l| can be found in the literature where the den- 
sity of states is sometimes called the spectral density. 
From its definition, it is clear that n(S) is positive. As- 
suming that the measure for the links is normalized to 1, 
the partition function at j3 = is 1. It can be shown [42j 
that, if the lattice has even number of sites in each di- 
rection, and if the gauge group contains — 1, then (3cos9 p 
goes into ~f3cos9 p by a change of variables 8i — > 9i + ir 
on a set of links such that for any plaquette, exactly one 
link of the set belongs to that plaquette. Using 

Z(-P) =c 2 ^Z(p) , (7) 
we find the duality 

n{2M p -S) = n(S) . (8) 

This implies the reflection symmetry 

/( S )=/(2-s). (9) 

Numerical values of f(s) have been obtained for dis- 
crete values of s between and 1. When s is close to 
or 2, f(s) diverges logarithmically and we can only reach 
values of s that are not too close to or 2. Consequently, 
the results cannot be used if |/3| is too large. Using the 
symmetry Eq. ^ and interpolation methods, a contin- 
uous function can be obtained in an interval [5, 2 — 5], 
where 6 is an appropriately small quantity. 

III. CALCULATION OF THE DENSITY OF 
STATES 

We performed Monte Carlo simulations in pure U(l) 
gauge theory using Biased Metropolis-Heatbath Updates 
[37l | . To cover a large range of couplings p 6 [0, 9] we used 
the multicanonical (MUCA) algorithm (3f| with Wang- 
Landau recursion [36| for the multicanonical weights. 
The software we used is described in Ref. [H|. 

We generated large statistics on symmetric lattices 
with volumes 4 4 , 6 4 and 8 4 . After the initial recursion 
we performed three MUCA runs on 4 4 , and two runs on 
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FIG. 1. f(s) on a 4 4 (top) , 6 4 (middle) and 8 4 (bottom) 
lattices; the errors have been multiplied by M p \ for readability, 
arbitrary constants have been added to separate the data sets 
and only one of every 40 points are displayed. 
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FIG. 2. Difference between the average plaquette calculated 
with the density of states and directly, for a 8 4 lattice. 



6 4 and 8 4 . The first MUCA run on 4 4 was regarded as 
exploratory and we did not include it in the final analy- 
sis. The weights for each next run were refined from the 
previous run. In total we used 20 independent streams 
for each lattice volume. In each stream we ran Wang- 
Landau recursion for the multicanonical weights before 
the production, therefore the weights differ between the 
streams, Wij(S), where S is the total action, i = 1, 20 
denotes different streams and j = 1,2 denotes MUCA 
runs. 

The quality of a MUCA run is indicated by the num- 
ber of tunneling events (i.e., how often during a run the 
system travels from the lowest energy to the highest and 
back). Also, to estimate how many statistically indepen- 
dent events we generated, we measured the integrated 
autocorrelation times. These data are summarized in 
Appendix Our statistics consists of N equ i sweeps for 
equilibration and N rpt = 64 x N equ i sweeps for measure- 
ments, where N equ i = 10 6 for 4 4 and 6 4 lattices and 
8 x 10 5 for 8 4 . 

For the error analysis we considered two MUCA runs 
in each stream as independent measurements. Thus, on 
each lattice we had 40 independent measurements in to- 
tal. For all quantities in the following the error bars are 
estimated from an uncorrclated average of these 40 mea- 
surements, weighted with the number of tunneling events 
in each corresponding run, since runs with more tunnel- 
ings sample the density of states better. The average 
results for f(s) are shown in Fig. [TJ 

To reweight an observable to the canonical ensemble 
we need to cancel the multicanonical weight Wij(S) and 
replace it with the Boltzmann factor exp(— f3S). For 
an observable O of interest, for instance, plaquette, we 
reweight the time series accumulated during a MUCA 
run ij to a given coupling (3: 



(0)ij(J3) = 



The final average is then given as 

^/yP' v-2 Ntunn ' 

where the number of tunnelings jV™ n " is given in Ap- 
pendix [A] 



IV. CONSISTENCY WITH EXISTING RESULTS 
AND EXPANSIONS 

A. Comparison with the average plaquette 

As a check of consistency, we compared the average 
plaquette at various /3, as obtained directly from the runs, 
Eq. (fTTj) . and calculated using the average density of 
states. As shown in Fig. [2j there is a good agreement 
within the estimated errors. 



B. Series for f(s) 

We compared the numerical results for /(s) with an- 
alytical results obtained using the weak and strong cou- 
pling expansions. The general methodology has been dis- 
cussed for SU (2) in [43| and remains applicable here. The 
basic ingredient is the saddle point equation at sq: 



/'(so) = P , 



(12) 



(10) 



which can be used to convert an expansion of / in powers 
of s (or (s — l) 2 ) into an expansion of sq in powers of 1//3 
(or (3, respectively). The coefficients of / can then be 
determined whenever the appropriate expansion of the 
average plaquette is available. In order to take the finite 
volume effects into account, we need to include at least 
the lowest order volume correction, namely 



P 



S o + (l/2A^)(/'"( So )/(/"( S o)) 2 ) + 0(l/AA p 2 ) . (13) 
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Using Eqs. (jT2|) and ([13)) together with an existing ex- 
pansion of P including l/J\f p effects up to a certain order, 
one can determine the coefficients of / up to the corre- 
sponding order. 

It should also be noted that at finite volume, there is 
a finite range of /? for which /' has a "Maxwellian kink" 
(discussed in Sec. IVIII) and three solutions of Eq. (fT2)) arc 
available rather than one. In this region, both expansions 
are expected to fail. We now proceed to discuss them 
separately. 



C. Strong Coupling 



U(l) 8 strong coupling 



order 2 
order 4 
order 6 
order 
order 10 
order 12 
num. error 




Following Rcf. |43[, we define 



g(y) = f(l + y)=Y,92r 



y 



2m 



(14) 



Using saddle point approximation and comparing, order 
by order, with the expansion of the average plaquette 
from the strong coupling expansion J44|, we can deter- 
mine the expansion of g(y). As in the case of SU(2), 
there are logarithmic singularities at y = ±1, which can 
be subtracted by defining 



h(y) = g(y) - A\n(l - y') 



2\ — 



E h 2m y 2m 



(15) 



The value of A comes from the weak coupling expan- 
sion and will be discussed in the next subsection (see 
Eq. PO]) ). The infinite volume results are summarized 
in Table |U The entries make clear that as the order 
increases, the effect of the logarithmic subtractions be- 
comes smaller. This indicates singularities closer to y = 

{8=1). 
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FIG. 3. Natural logarithm of the absolute value of the dif- 
ference between f(s) calculated on a 8 4 lattice and successive 
approximations obtained from the strong coupling expansion. 



Indeed, they are even hard to resolve for V = 4 4 . This 
can be traced to the fact 43| that even for this volume, 
the dependence on V would appear at order /3 8 from the 
contributions of strong coupling graphs called torelons 
[45j that wrap around the periodic volume in one di- 
rection. As translations in that direction do not gener- 
ate new graphs, such graphs have a suppression of order 
l/L compared to graphs with a trivial topology. Con- 
sequently, for order less than 8, the finite volume effects 
can be estimated by canceling the volume dependence in 
the two terms of the r.h.s. of Eq. ([T3|) . For instance at 
lowest order, we find that 



g 2 = -1 + 1/(810 



(16) 



Consequently the difference A/(s) of /(s) for two differ- 
ent volume V\ and V 2 near y = (s = 1) is 



A/( s )c(l/8)(l/U 1 -l/U 2 )( S -l) 5 



(17) 



Even for V = 4 4 , this difference is smaller than the error 
bars in the region s ~ 1. In order to reduce the noise, we 
have averaged the data in bins of ten data points. The 
results arc displayed in Fig. U which shows that the data 
and the analytical result in Eq. (| 1 T[) arc compatible. 



D. Weak coupling 



TABLE I. U(l) strong coupling expansion coefficients a2 m of 
P (rescaled from Ref. 13), and of f(s) defined in the text. 

The improvement of the approximation with successive 
order is shown in Fig. [3J The graph shows that for 
s ^ 0.5, successive orders provide better approximations 
up to the point where the numerical accuracy is reached. 
Such a range corresponds to a convergent region j3 ^5 0.9 
in the /3-plane. 

It should be noted that for the strong coupling expan- 
sion, the finite volume effects are negligible for V = 8 4 . 



A similar approach can be followed in the weak cou- 
pling limit. At small s, the logarithmic singularity dom- 
inates and we assume that 



f(s)=Aln(s)+J2fmS r ' 



(18) 



m=0 



The unknown coefficients can be determined from the 
weak coupling expansion of the average plaquette 



p ~ Yl 



(19) 
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U(l) 8 weak coupling 
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FIG. 4. Difference between f(s) on 4 4 and 8 4 near s =1 FIG. 5. Difference between f(s) calculated on a 8 4 lattice and 

(boxes). The circles are obtained by averaging over bins of successive approximations obtained from the weak coupling 

size 10. The solid line is Eq. (|17[) . The part with s > 1 can expansion, 
be obtained by symmetry and is not shown. 



The volume dependent coefficients b m have been calcu- 
lated up to order 4 in Ref. [46| . The two lowest orders 
of the expansion can be performed exactly and yield: 

A = l/4-5/(12V) , (20) 
h = (!/8)(l - 1/V) . (21) 

The higher orders involve numerical loop calculations. 
The results of the expansion as well as the volume cor- 
rections for the 4 4 and 6 4 lattices are shown in Table HT1 





V = oo 


V = 6 4 


V = 4 4 


h 
b 2 
bz 


i 

4 
1 


1295 
5184 
2171747375 


255 
1024 
65025 


32 

0.01311 


208971104256 

0.01309 


2097152 

0.01296 


bi 


0.00752 


0.00749 


0.00739 


A 

h 

h 


i 

4 
1 
8 

0.07363 


3883 
15552 

1295 
10368 

0.07359 


763 
3072 

255 
2048 

0.07314 


h 


0.07638 


0.07605 


0.07515 



TABLE II. The weak coupling expansion for V = oo, 4 4 and 
6 4 . The upper half is the list of the expansion coefficients of 
the average plaquette P with respect to 1//3 [46||. The lower 
half is the corresponding list of expansion coefficients of f(x). 

The difference between numerical and analytical re- 
sults is shown in Fig. [5] for V — 8 4 . The graph makes it 
clear that the quality of the approximation increases as 
three successive corrections to the leading logarithm are 
added. The third correction is good enough to reproduce 
the data within the numerical accuracy for s < 0.1. This 
order is not sufficient to identify a "non-perturbative en- 
velope" defined in Ref. [47| and observed for SU(2) in 
Ref. [13. 



U(l) 4 4 minus 8 4 Weak Coupling 
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FIG. 6. Difference between /(s) on 4 4 and 8 4 near s = with 
average over bins of size 10. The solid line is the expansion 
described in the text. 



The difference between V — 4 4 and V — 8 4 is shown 
in Fig. [5] Note that for the smallest volume (V = 4 4 ), 
the resolution in s used during the multicanonical sim- 
ulation is coarser than the 1,000 bins used to represent 
f(s). Consequently, some small "staircase" structure ap- 
pears near where /(s) changes rapidly. For this reason, 
we have averaged f(s) over bins of size 10 and Fig. [5] 
shows a good agreement with the analytical expansion 
that includes the logarithmic singularity and a linear 
term. Higher order corrections are significantly smaller 
than the errors bars. There is an arbitrary constant in 
the expansion of f(x) which cannot be determined by 
the saddle point equation. For the numerical data, such 
a constant may differ for two different volumes and needs 
to be subtracted. 
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U(l) 4 4 minus 8 4 Crossover 
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FIG. 7. Difference between f(s) on 4 4 and 8 4 near s — 0.5. 
For readibilty, we only show every 5 points (no binning). The 
solid line is the fit given in Eq. (|22[) . 



U(l) 6 4 




.34 0.36 0.38 0.4 0.42 
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FIG. 8. f(s) - j3s for /3 = 1.00175, 100177 and 1.00179 on a 
6 4 lattice. The horizontal lines have been drawn to emphasize 
small height asymmetries. 



V. VOLUME DEPENDENCE IN THE 
CROSSOVER REGION 

A. Empirical parametrization 

The difference of f(s) for 4 4 and 8 4 resembles the ef- 
fective potential for the central Coulomb potential with 
a leading singularity near s = 0.35 and a 1/s behavior at 
larger s. Using in addition a constant that has no partic- 
ular meaning as long as we don't normalize the density of 
states and a 1/s 2 correction, we performed a 4-parameter 
fit with the 311 bins corresponding to 0.39 < s < 0.7. 
The numerical result is: 

f(s) = - 0.00112063 + 4.82641 x 10" 6 /(-0-35 + sf 
- 0.000680501/s 2 + 0.00172882/s. (22) 

As shown in Fig. [3 it fits the data reasonably well. 

B. Volume dependence of the double peak 

The plaquctte distributions for the volumes considered 
here have a double peak structure for j3 near 1. At finite 
volume, it is easy to locate the value of j3, denoted Ps 
hereafter, where the two peaks of f(s) — Pss have equal 
height. Other pseudocritical (3 have been defined in the 
literature [3(1 |4§ - |50| . The accuracy of the determina- 
tion of (3s depends on the smoothness of the distribution 
and the size of the error bars. In Fig. El we show that 
f(s) - j3s is slightly tilted to the left for (3 = 1.00175 
and to the right for (3 = 1.00179. Given the smoothness 
of the distribution, we conclude that (3s = 1-00177(2). 
With the same graphs, we can also determine approxi- 
mate values of the two maxima si and S2- The numerical 
results for the three volume considered are provided in 
Table HH 

The density of states can be used to calculate the pla- 
quette probability distribution at /3s- The results are 



L 


Ps 


Sl 


S2 


4 


0.9793(1) 


0.370(5) 


0.445(5) 


6 


1.00177(2) 


0.353(2) 


0.411(2) 


8 


1.00734(1) 


0.349(1) 


0.395(1) 



TABLE III. /3s, si and s 2 defined in the text for L = 4, 6 
and 8. 



shown in Fig. [9] where the normalization has been cho- 
sen in such a way that the integral under the curve is 
approximately one. Fig. |9] makes it clear that the dip 
between the peaks deepens and the peak separation de- 
creases as the volume increases. This will be discussed 
in more detail in Sec. IVIII 



U(l) plaquette distr. at p s 
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FIG. 9. Plaquette distribution for U(l) at Ps for L= 4, 6 
and 8. 
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VI. FISHER'S ZEROS 



A. Approximate zeros from reweightings 

Approximate values of Z at fixed /3 can be obtained by 
using the Ricmann sum approximations of Eq. (j4|) : 



U(l) 4 ln|8z/Z 



(23) 



We can now study how the error Sf on f(x) can affect 
our estimates of Z. The relevant quantity is Af p 5f is 
included in Fig. [T] For the three volumes, N p 8f is of 
the order of a few percents and linearization is justified. 
Small scale fluctuations of the same order are visible in 
the distributions of the independent streams. 

As we are interested in locating Fisher's zeros, it is 
clear that the errors have a potentially important effect 
near an approximate zero. The best we can do is to 
identify regions where \Z\ is significantly larger than \SZ\ 
so that we can confidently say that there are no zeros in 
these regions. If we use the linear estimate 



we have the inequality 



\SZ(p)\<AsJ2K\8f(s)\e 



X p {f(s)-Re{ts) 



(24) 



(25) 



but in general the bound is not sharp because the sign 
of 5f(s) can vary rapidly. We have estimated \5Z\ by 
taking the difference between Z calculated with the av- 
eraged / and Z calculated with the stream with the most 
tunnclings. The results are shown in Fig. 1101 A mesh of 
0.00025 in /3 is used which is larger than the typical fluc- 
tuations in /. The light (toward yellow on-line) regions 
represent the areas where we cannot exclude zeros. The 
dark (toward blue on-line) regions represent the areas 
where zeros are very unlikely. A small light region is an 
indication for the existence of a zero while a broad light 
region indicates that the errors dominate. The second 
possibility typically appears at large imaginary (3 where 
due to rapid oscillations of the integrand, cancellations 
occur making the final results more sensitive to the errors 
on f(s). In view of this remark, Fig. [T0l suggests that 
reweighting methods allow to estimate the locations of 
the two lowest zeros for L = 4 and three lowest zeros for 
L = 6. 

We have also calculated KcZ and ImZ from Eq. (f23|) 
using the average /. Their respective zeros are shown in 
Fig. QT] The complex zeros appear at the intersections 
of the two sets of curves defined by ReZ=0 and ImZ=0 
respectively. This happens in a way which is consistent 
with Fig. 1101 Error bars can be estimated by comparing 
the intersections for the 20 streams. The results are given 
in Tables HVl and IV1 
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FIG. 10. \n\SZ/Z\ for 4 4 and 6 4 lattices (color online). 



L 


first 


second 


third 


4 


0.9791(1) 


0.9780(4) 


not stable 


6 


1.00180(5) 


1.0007(1) 


0.9993(5) 


8 


1.00744(2) 


1.0068(2) 


1.0061(4) 



TABLE IV. Real part of the first three zeros. 



L 


first 


second 


third 


4 


0.0259(1) 


0.057(1) 


not stable 


6 


0.00758(2) 


0.018(1) 


0.027(2) 


8 


0.00306(2) 


0.008(1) 


0.012(1) 



TABLE V. Imaginary part of the first three zeros. 
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FIG. 11. Zeros of the Re (+, blue on-line) and Im (x, red 
on-line) part of Z for (7(1) using the density of states for 4 4 
and 6 4 lattices. 



B. Chebyshev interpolations 

The original grids of the density of states are sometimes 
not sufficient for precise numerical integrations which is 
how we define our partition function. It is especially 
true when the imaginary component of f3 is large and, 
as a consequence, the partition function oscillates more 
frequently than the original grids can resolve. It is con- 
venient to apply the Chebyshev interpolation which pro- 
vides arbitrary integrating step sizes for designed integral 
precision. For the Chebyshev interpolation of numerical 
data, the determination of the coefficients by the Least 
Square Fit method is more efficient and robust than by 
discrete or integration methods. In this paper, we will 
primarily follow this approach. 

A range of interest [a, b] can be mapped to [—1, 1] in 
which we express the target function by 



N r . 



We then minimize the distance of the function to a data 
set or multiple data sets, which will uniquely determine 
the coefficients c„ of linear equations. 

We shall keep in mind that, like other polynomial ap- 
proximations, Chebyshev interpolations may introduce 
artifacts such as fake zeros. We want to make sure that 
the true zeros are distinct from the fake ones. Special at- 
tention should be paid to the range of approximation. 
In practice, we often use a small range to emphasize 
the numerical signal from a certain region. The aver- 
age plaquette (x) is related to the coupling f3, through 
(x) fj = -d\nZ{p)/dp/N p . This is not valid if (x) goes 
beyond the range of the approximation (an ellipse in the 
complex plane, see below). Reducing the range of inter- 
polation may introduce fake zeros with large Im/? in the 
f3 complex plane. However the lowest zeros are usually 
not affected. Care should also be paid on the orders of 
the Chebyshev approximation. True zeros should be in- 
dependent of the order of polynomials. In the following, 
we always use various orders of Chebyshev interpolations 
and make sure that the zeros arc free of these artifacts. 



C. Ellipse of convergence 



The definition 



T n {z) = cos(n arcos(z)) 



(27) 



shows that the expansion in Chebyshev polynomials is a 
Fourier expansion for the variable arccos(x). If |c„+i/c n | 
from Eq. ([26]) reaches a limit C, then the expansion con- 
verges for |T„(z)| < C~ n . To work on the complex plane, 
the following relation is helpful: T n (z) = (uj n + co~ n )/2, 
when z is expressed as z = (lo + lu~ 1 )/2. The conver- 
gence of a Chebyshev series is then analyzed through the 
variable u>. It can be shown [5ll | that the region of con- 
vergence on the w-plane is a ring confined by a pair of 
concentric circles and the region is mapped into an area 
bounded by an ellipse on the z-planc. 

The continuation of the Chebyshev expansion to the 
complex plane is limited by the ellipse. Fortunately, in 
the case of U(l), the lowest complex zeros are typically 
very close to the real axis and these zeros are well inside 
the ellipse of convergence. 



D. Locating zeros with the residue theorem 

There is a general algorithm to find the zeros of an 
analytic function by using Cauchy's Integral Theorem 
[52^ . For simplicity, we will only consider the special case 
when all the zeros are of order 1 which apparently applies 
to our problem. Suppose that an analytic function Z(JS) 
has K zeros enclosed by a closed contour C, then 



f(y) = z^c n T n {y) 



(26) 



^-.j{hxZ)'p n dp = Y,{PiT, n = 0,1,2, 



(28) 
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L 


Re/3 




dc 


Im/3 


(J s 


Cc 


4 


0.9791235 
0.9777314 
0.9752954 


3.6e-5 
3.5e-4 
l.le-3 


5.3e-8 
7.1e-6 
2.9e-4 


0.0260065 
0.0572764 
0.0831705 


3.7e-5 
1.4e-4 
1.3e-3 


3.9e-9 
3.3e-6 
3.2e-4 


6 


1.0017969 
1.0007433 
0.9988964 


1.7e-5 
6.0e-5 
1.4e-4 


1.7e-6 
2.3e-5 
2.7e-4 


0.0075821 
0.0182044 
0.0271866 


8.7e-6 
2.8e-5 
4.5e-4 


1.4e-6 
4.0e-6 
1.5e-4 




1.0074380 
1.0068296 
1.0060410 


l.le-5 
2.3e-5 
l.le-4 


7.7e-8 
2.1e-6 
1.2e-5 


0.0030653 
0.0077673 
0.0115079 


3.6e-6 
2.4e-5 
1.0e-4 


6.8e-8 
3.3e-7 
8.5e-5 



U(l) 



TABLE VI. The lowest three zeros in the volumes 4 4 , 6 4 and 
8 4 . Columns 2-4 are the real parts of the zeros, the estimate 
error a s from different streams of Monte Carlo runs and the 
error a c due to the orders of Chebyshev interpolation (we used 
three different orders 40, 44 and 50 for all three volumes). 
Columns 5-7 are similar quantities for the imaginary parts. 



where Pi are all the zeros in contour C. When n = 0, the 
summation on the right hand side is just the number of 
zeros. 

The partition function we are considering is an analytic 
function, since it is just a sum of analytic functions. We 
scan the complex plane with rectangular contours which 
enclose two or less zeros. We monitor the n = integral 
which should give the number of zeros very close to an 
integer and a very small imaginary part. The method 
turns out to be quite robust and reliable. 



E. Zero structure near the real axis 

The lowest zeros from three volumes are given in Ta- 
ble I VII and shown in Fig. [T^l The error bars take into 
account both the Monte Carlo statistical error and the 
(much smaller) Chebyshev interpolation error. The three 
lines are the linear fits for the first, second and third low- 
est zeros. They intersect the real axis approximately at 
the same point f3 = 1.01134(1). Fig. [T^Jalso shows that 
Ps and the real part of the zeros are highly correlated. 

The good look of the linear fits is deceptive as they 
have a rather large x 2 an d a small goodness of fit Q (see 
p. Ill of [53[) which can be explained by the small er- 
rors bars. Another potentially deceptive result is that the 
imaginary part of the lowest zero decreases like L~ 3 08 . 
If this result is indicative of what happens at larger vol- 
ume, this would be interpreted as signaling a second order 
phase transition with v ~ 1/3.08 ~ 0.325. Larger lattices 
are needed as will be discussed in Sec. IVII1 



o . 1 



.08 



.06 



0.04 



. 02 





ii 

1 Qdst) = 


ill 

J 075 


: 


Q(2d) = 


.068 - - 




Q(3rd) < 


.001 


- 








\ 


6 4 ™ ' 






8 4 " - 


- 


\\ 


- 


- 


\ 


- 






















,1,,, 


1 , , \, ^ 


96 


. 98 


1 1 




Re(p) 





FIG. 12. The lowest zeros from the volumes 4 , 6 and 
8 4 (from left to right). Linear fits for the first, second and 
third zeros (bottom to top) and their goodness of fit Q. The 
diamonds on the real axis are the double-peak /3's from Table 

imi 



L 


Sa 


Sb 


Pa 




4 


0.274 


0.488 


1.125 


0.945 


6 


0.284 


0.436 


1.1 


0.985 


<s 


0.295 


0.408 


1.075 


1 



TABLE VII. Values of s a and st and the corresponding values 
of/3 (P(p) = s). 



zeros near j3 = 1 , the density of states in a finite range is 
sufficient. This information becomes very important at 
higher volumes where the calculation is more expensive. 
In Table IVII1 we provide the values s a and Sb outside of 
which the knowledge of f(s) has effects smaller than the 
error bars on the zeros. 



VII. TOWARD LARGER VOLUME 
CALCULATIONS AND SCALING 

In this section, we explore ways of discriminating be- 
tween first and second order transitions. For this pur- 
pose, we include data at larger volumes (L=10, 12, 14 
and 20) from Ref. [56| with much lower statistics, namely 
one stream with two MUCA runs. 



Zeros 



F. Dependence on the range of integration 

In the previous calculations, the tails of integration 
play a marginal role. If we are interested only in the 



In Sec. I VII we explained that for the data for L = 4, 6 
and 8, the imaginary part of the lowest zeros scales like 
£-3.08 j|- j s p 0SS ible that as the volume increases, the 
approach of the real axis "rolls" toward the L~ A scaling 



10 



L 


Re/3 


ARe/2 


Im/3 


Aim/ '2 


10 


1.00947 


2e-5 


0.001478 


2e-6 


12 


1.01027 


2e-5 


0.000795 


2e-6 


14 


1.01064 


2e-5 


0.000449 


8e-6 


20 


1.01101 


le-5 


0.000119 


le-6 



TABLE VIII. Higher volume zeros. Columns 2 and 4 are the 
averages over the two MUCA runs. Columns 3 and 5 are one 
half of the differences (not the estimated error, see text). 



expected for a first order transition [30[ • We now discuss 
the scaling of the zeros using the lower statistics data for 
the larger volumes given in Tabic IVIII1 

It is questionable that two MUCA runs could lead to 
a reliable estimate of the errors. An error bar from just 
two independent measurements fluctuates strongly and 
reaches a 95% confidence range only at about 14 (instead 
of 2) error bars (see p. 78 of [531 )■ We decided therefore 
to smoothen the error bars by assuming that the real rel- 
ative error is the same for all four of our large lattices. 
Averaging these relative errors and multiplying the by 
three, the approximate 95% confidence range of four in- 
dependent data, gives an error bar of 1.69%, which is 
then given in the fourth column for all the large data 
of of Tabic IIXI Not to overweight the far more accurate 
small lattice against the large lattice data in the sub- 
sequent fits, they are also used with a relative error of 
1.69% and thus listed in Table llXl We want to empha- 
size that this procedure has been designed to understand 
how different fits allow to discriminate between first and 
second order rather than to extract accurate values for 
the fitted parameters. 



L 


First Run 


Second Run 


Combined 


4 






0.02691 (44) 


6 






0.00758 (13) 


8 






0.003065 (52) 


10 


0.0014756 


0.0014797 


0.001478 (25) 


12 


0.0007927 


0.0007969 


0.000795 (14) 


14 


0.00045747 


0.00044157 


0.0004495 (76) 


20 


0.000011882 


0.000011901 


0.00001189 (21) 



TABLE IX. y — Imz from two independent runs on L > 
10 lattices and their combination as explained in the text 
together with reduced accuracy values from L < 8 lattices. 



For these data we performed 3-parameter fits, which 
are listed in the following together with their goodness 



0.01 



0.001 



0.0001 



Fit L=4,6,8 
Fit L=4 to 20 



Q=0.39 



Q=0.43 
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FIG. 13. Fits of lmz(L) on a log-log scale. 



of fit Q. 



ai_ 
L 4 



1 



y = a x L a * ( 



«2 »3 

L L 2 

, "3 



Q = 0.43, 



j-) , Q = 6.2 x 10~ 4 



y 



L 4 



(l + a 2 L a3 ), Q = 2.8x10" 



(29) 
(30) 
(31) 



The first fit shows that L~ 4 behavior is consistent with 
all the data put together. The other two fits are in dis- 
agreement with the data. 

Using only the data from the L = 4, 6, 8 lattices with 
the modified error bars given in Table llXl the 2-parameter 
fit 



V = ai L° 



Q = 0.39 



(32) 



is also in agreement with the data and gives the exponent 
a 2 = -3.082 (35) instead of -4. The fits ^ and 
are shown in Fig. 1131 
However, the fit 



!J 



£3.08 



1 



a 2 

T 



03 

L 2 



(33) 



with the seven data points leads to Q < 10~ 8 . In addi- 
tion, if we perform a four parameter fit as in Eqs. (|29[) 
and (|33p but with the leading exponent fitted, we ob- 
tain 4.121(74) for this exponent with Q=0.72. These 
results seem to favor the first order possibility However, 
they should be checked with higher statistics data for the 
larger volumes. 



B. Features of f(s) 

In the infinite volume limit, the width of the double 
peak distribution goes to a nonzero limit (latent heat) 
for a first order phase transition. For a second order 
transition, this width should go to zero as an inverse 
power of L. These two possibilities are tested by plotting 
the width versus 1/L or in a log-log scale in Fig. Q31 
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FIG. 14. Width of / - p s s as a function of L. 



The difference between the local minimum and the lo- 
cal maxima of f — Pss (by definition of /3s, the two local 
maxima have the same height) should decay like C/L for 
a first order transition with C proportional to the inter- 
face tension. For a second order transition, this difference 
should go to zero as an inverse power of L. The data is 
shown versus 1/L and on a log- log scale together with 
simple fits (made without L = 4) in Fig. [15] In the first 
fit, the 1/L 2 corrections are clearly important and it is 
not surprising that the power in the second fit is between 
-1 and -2. 

Using arguments from Ref . [lij , leads to the conclusion 
that for a second order phase transition, the width should 
scale like L - ' 1- ")/^ , while the depth should scale like 
L~ D . If we use D = 4, v ~ 0.325 and the hyperscaling 
relation a = 2 — Dv ~ 0.7, we obtain (1 — oi)/v ~ 0.92 
which is not too far off for the width but very far off for 
the depth. 



VIII. CONCLUSION 

Using multicanonical methods we have calculated the 
density of states for pure U(l) lattice gauge theory with 
high precision on small 4 4 , 6 4 and 8 4 lattices and with 



Depth of dip in f(s)-p s * s 
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1/L 
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FIG. 15. Difference between the local minimum and the local 
maxima of / — ftss as a function of L. 

moderate precision on larger 10 4 , 12 4 , 14 4 and 20 4 lat- 
tices. From these data we were able to locate low-lying 
Fisher's zeros by Chebyshev interpolations and residue 
theorem methods. On the small lattices the scaling prop- 
erties of the zeros are consistent with a second order 
phase transition, while from the larger lattices there is 
some indication that this turns around and becomes con- 
sistent with a first order transition. 

Although £7(1) lattice gauge theory was already in- 
troduced in the pioneering paper by Wilson (55|, it still 
resists to reveal clearly the true nature of its transition 
from the confinement to the Coulomb phase. Like other 
physical quantities, e.g., Polyakov loop susceptibilities, 
Fisher's zeros appear to need rather large lattices to dis- 
play their asymptotic scaling properties. As modern su- 
percomputers allow parallel processing on an unprece- 
dented scale, the solution may finally become achieved 
by brute force calculations on very large lattices. 
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U(l) 6 Im. Part of Zero 




0.00756 0.00758 0.0076 

stream 1 



0(1) 6 4 Real Part of Zero 




1.00175 1.0018 1.00185 

stream 1 





MUCA 1 


MUCA 2 


MUCA 3 


# 






Tint 


j^tunn 


Tint 


1 


1,512 


3,886 


99(2) 


4,697 


110(1) 


2 


1,406 


4,018 


109(1) 


2,833 


720(35) 


3 


1,974 


4,723 


114(5) 


3,492 


540(15) 


4 


1,552 


4,537 


216(5) 


2,684 


750(18) 


5 


774 


5,038 


175(12) 


4,920 


367(29) 


6 


963 


4,769 


109(1) 


5,682 


164(4) 


7 


196 


397 


1383(124) 


3,162 


679(34) 


8 


4,089 


3,875 


101(5) 


4,547 


116(5) 


9 


2,344 


4,214 


108(7) 


4,599 


266(6) 


10 


1,652 


4,582 


185(15) 


3,484 


625(48) 


11 


1,622 


4,179 


111(7) 


5,030 


255(12) 


12 


1,722 


4,281 


126(3) 


3,985 


674(101) 


13 


3,406 


4,081 


98(1) 


4,994 


146(3) 


14 


1,271 


4,127 


104(6) 


5,257 


135(6) 


15 


488 


4,610 


255(9) 


3,776 


524(14) 


16 


2,351 


4,394 


108(3) 


4,167 


338(12) 


17 


788 


4,785 


123(4) 


4,598 


356(10) 


18 


735 


4,680 


134(4) 


4,661 


364(16) 


19 


845 


4,450 


200(4) 


3,675 


537(14) 


20 


2,697 


3,526 


93(1) 


4,123 


104(7) 



TABLE X. MUCA data for 4 4 . 



FIG. 16. Imaginary (top) and Real (bottom) part of the com- 
plex zeros of MUCA run 1 versus MUCA run 2 for the 20 
streams on a 6 4 lattice. 
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Appendix A: Multicanonical data 

The number of tunnelings and the integrated autocor- 
relation times for 4 4 , 6 4 and 8 4 lattices is given in Ta- 
bles IXllXIII Fig. [16] illustrates the fluctuations among 
streams and the correlations among MUCA runs for the 
zeros on a 6 4 lattice. In table IXlLTl we summarize the 
parameters of simulations on 10 4 -20 4 lattices. 



Appendix B: Numerical data used in Sec. IVIII 





MUCA 1 


MUCA 2 


# 




Tint 


jn^tunn 


Tint 


1 


1,351 


702(48) 


950 


559(27) 


2 


708 


466(12) 


897 


480(11) 


3 


367 


422(12) 


947 


502(17) 


4 


454 


474(30) 


971 


561(31) 


5 


580 


484(37) 


894 


565(47) 


6 


682 


481(21) 


911 


511(23) 


7 


523 


423(14) 


909 


557(43) 


8 


765 


512(28) 


903 


532(33) 


9 


696 


510(27) 


921 


485(12) 


10 


652 


469(32) 


935 


569(53) 


11 


513 


544(46) 


976 


523(33) 


12 


378 


396(18) 


976 


536(27) 


13 


867 


559(17) 


961 


495(10) 


14 


661 


496(13) 


932 


509(28) 


15 


545 


542(50) 


896 


554(38) 


16 


615 


497(14) 


920 


496(13) 


17 


475 


438(16) 


979 


543(46) 


18 


822 


464(10) 


903 


486(11) 


19 


878 


570(35) 


892 


508(30) 


20 


578 


588(62) 


949 


572(28) 



In order to calculate the quantities used to make the 
graphs of Sec. IVIII we have first fitted f(s) in a re- 
gion slightly wider than the location of the two peaks 



TABLE XL MUCA data for 6 4 . 



13 



L 


Ps 


Sl 


S2 


A(/-/M 


sl 


SR 


Pl 


Pr 


Pr-Pl 


4 


0.979327 


0.369215 


0.442409 


0.0000597527 


0.384952 


0.425551 


0.976759 


0.981792 


0.00503262 


6 


1.00176 


0.352413 


0.410538 


0.0000533197 


0.364762 


0.395541 


0.998768 


1.00446 


0.00569344 


8 


1.00738 


0.348089 


0.394932 


0.0000337064 


0.358049 


0.382876 


1.00504 


1.00951 


0.00447297 


10 


1.00942 


0.345567 


0.386254 


0.000023194 


0.353401 


0.37509 


1.00747 


1.01103 


0.00356233 


12 


1.01022 


0.345633 


0.380806 


0.0000164495 


0.352382 


0.370685 


1.0086 


1.01153 


0.00293064 


14 


1.01058 


0.345352 


0.377624 


0.0000114541 


0.351326 


0.368136 


1.00934 


1.01155 


0.00221252 


20 


1.01097 


0.345337 


0.374469 


7.1542782 10 -6 


0.349468 


0.36223 


1.00988 


1.01157 


0.0016911 



TABLE XIV. Numerical results described in the text. 





MUCA 1 


MUCA 2 


# 


j\jtunn 


Tint 




Tint 


1 


145 


1,756(99) 


252 


1,782(124) 


2 


79 


1,298(69) 


240 


1,798(78) 


3 


75 


1,475(144) 


259 


2,486(368) 


4 


121 


1,291(141) 


216 


2,159(235) 


5 


150 


2,420(364) 


216 


1,660(154) 


6 


86 


1,097(40) 


256 


1,528(55) 


7 


74 


1,103(51) 


255 


1,621(105) 


8 


132 


1,419(48) 


254 


1,568(98) 


9 


187 


3,089(438) 


197 


3,074(397) 


10 


98 


1,370(94) 


255 


1,706(104) 


11 


142 


1,389(47) 


208 


1,700(218) 


12 


165 


1,804(343) 


254 


1,800(256) 


13 


93 


1,265(111) 


270 


1,912(160) 


14 


212 


2,012(114) 


231 


1,855(254) 


15 


137 


1,581(181) 


249 


1,855(170) 


16 


159 


1,904(235) 


211 


1,674(224) 


17 


269 


1,773(76) 


213 


2,169(256) 


18 


206 


1,773(176) 


234 


1,503(52) 


19 


214 


1,756(76) 


212 


1,954(265) 


20 


96 


1,680(221) 


227 


1,697(101) 



volume 


sweeps 


ftmin, 


ftmax 


MUCA1 


MUCA2 


10 4 


32 x 96, 000 


0.980 


1.030 


103 


133 


12 4 


32 x 112, 000 


0.990 


1.030 


75 


82 


14 4 


32 x 128, 000 


1.000 


1.020 


57 


51 


20 4 


64 x 100, 000 


1.010 


1.012 


155 


210 



TABLE XIII. MUCA data for volumes where simulations were 
performed in narrow [/3 m ; n , fimax] range. Last two columns 



summarize the number of tunneling events. 10 4 
from 15611. 



12" 



TABLE XII. MUCA data for 



of /(s) — fls with the first 12 Chebyshev polynomials. 
Higher order polynomials tend to pick up the noise and 
provide results which are less stable. Using this polyno- 
mial approximation, we calculated the two roots of f"(s) 
in the interval considered. We call them (Left) and s ^ 
(Right). They are the local extrema of f'(s). The cor- 
responding P (obtained from the saddle point Eq. (fT2]) 
are denoted j3 L and (3r. For /?£</?< (3r, the saddle 
point Eq. (fl"2f has three solutions instead of one (the 
"Maxwell kink"), Ps corresponds to the case where the 
area of the kink below and above are equal. The location 
of the maxima of f{s) — ftss are called si and S2 as in 
Sec. El A(/ — fiss) denotes the difference between the 
local minimum and the local maxima of / — /3s s. The 
numerical results are provided in Table IXIV1 
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